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ABSTRACT 

A Smooth  Particle  Hydrodynamics  (SPH)  method  has  been  developed  to  provide  a transient  structural  analysis 
of  fan  blades  during  bird  strikes  [I  J,  [3J 

The  bird  is  modeled  by  a set  of  lagrangian  particles  and  a specific  state  of  equation  evaluates  the  pressure  in 
the  projectile.  The  SPH  method  has  been  coupled  with  the  Finite  Elements  method  by  a contact  algorithm  and 
implemented  in  the  fast  dynamic  code  EUROPLEXUS  [2],  [4]. 

Unlike  the  classical  Lagrangian  formulations,  this  model  avoids  the  severe  distortions  of  projectile  mesh 
during  impact  phase.  An  uniform  bird  slicing  by  the  rotating  blades  is  not  assumed  prior  to  the  calculation. 
The  method  has  been  applied  to  compute  the  response  and  the  residual  deformations  of  a bladed  disk. 

1.  INTRODUCTION 

Bird  strike  in  a turbofan  engine  is  a phenomenon  which  can  occur  at  each  flight’s  stage  and  principally  during 
take  off  and  landing.  These  impacts  generate  severe  loads  on  the  fan  stage  and  sometimes  the  strains 
undergone  by  the  structure  can  lead  to  failure. 

Prior  to  the  final  certification  test,  some  expensive  spin  tests  are  usually  performed  to  validate  the  initial 
design  choices.  A failure  during  a bird  ingestion  will  be  very  expensive  because  of  the  necessary  re-design  of 
the  blade  (costs,  delays).  To  reduce  the  number  of  these  tests  and  to  avoid  any  delay  during  the  certification 
tests  the  manufacturers  perform  analyses  with  non-linear  transient  structural  codes. 

The  task  lies  then  in  the  development  of  a convenient  and  suitable  industrial  strategy  to  demonstrate  the 
reliability  of  a Finite  Element  (FE)  modelling  to  predict  the  impact  loads  and  damage. 

The  classical  Finite  Elements  method  do  not  allow  to  treat  completely  this  problem.  The  Lagrangian 
formulation  is  too  tricky  to  handle  for  the  bird  because  of  large  deformations.  The  coupled  Euler-Lagrange 
method  creates  a significant  numerical  dissipation  despite  difficult  remeshing  in  the  case  of  rotating  structures. 
Moreover,  CPU  time  widely  increases  for  these  last  calculations  because  a large  area  has  to  be  meshed  to 
anticipate  fluid  deformations  and  blades  rotation.  So,  all  these  methods  present  major  drawbacks  to  model  bird 
the  slicing  by  a rotating  structure. 

To  avoid  all  this  problems,  SNECMA  calculates  bird  strike  with  a SPH  method.  This  formulation  has  been 
implemented  in  the  fast  dynamic  analysis  code  EUROPLEXUS  developed  jointly  by  the  French 
“Commissariat  a TEnergie  Atomique”  and  by  the  JRC  of  the  European  Community  at  ISPRA. 

2.  THE  SPH  METHOD 

The  SPH  [5]  method  has  been  defined  by  astrophysicists  during  the  70’s.  But  at  the  beginning  of  the  90’s, 
some  mechanical  scientists  have  adapted  this  method  to  treat  impact  at  very  high  velocities  [6J,  [7J,  [8J. 

The  SPH  method  is  well-suited  to  simulate  problems  that  present  mesh  distortions  and  large  displacements. 
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2.1.  Mathematical  principles 

The  discrete  field  / .scalar  or  not,  defined  in  a domain  Q can  be  approached  and  regularised  by  a field  family 
(/)^ , defined  by  : 

{f)^  = \j{y)'W{x-y,h)-dy  (1) 

with  W Kernel  or  weight  function 

h smoothing  parameter  which  determine  the  volume  Q,  around  a:  where  W takes  its  greatest  values. 
This  parameter  h defines  the  regularity  of  (/)^ 

The  field /can  be  the  velocity,  the  density,  the  pressure  ... 

The  volume  Qi  is  smaller  than  the  domain  Q and  the  kernel  W generally  verify  the  following  properties  : 

1.  W{x-y,h)>Q  m QiW\\h  Qi  c Q 

2.  W(x-y,h)  = 0 VygQi 

3.  ^W{x-y,h)-dy  = \ (consistent  condition) 

4.  W(s,h)  is  monotonous  and  decreases  with  s = ||x-  y|| 

5.  W(s,h)  ^ S(s)  with  where  <?(.?)  is  the  Dirac  distribution 


The  most  common  kernels  W are  the  Gauss,  cubic  or  quadratic  fonctions. 

If  the  values  of  the  field /are  known  for  a finite  number  N of  points  jc„  the  integral  (1)  can  be  approximated  by 
a directe  sum  : 


(=1 


/(w) 

p(^) 


■W{x  — x.,h) 


(2) 


ftl 

The  mass  m,  is  associated  to  the  particles  x.  and  the  ratio  — — represents  the  volume  of  the  domain  Qi 

p(l) 

connected  to  Xi. 

The  equation  (2)  looks  like  form  functions  in  the  finite  element  formulation. 

The  approximation  of  the  field/can  be  easily  applied  to  the  fundamental  equations  of  Mechanics.  At  the  node 
number  k : 

■ continuity  equation  : 

^^^^^=p(x,J)■y—P^■(v{x,J)-v(x.,t))■V^W(\X|,  -x.\,h)  (3) 

at  M P(x,,t) 


■ momentum  equation  : 
dt 


ip)(x,,t)  {p)(x„t) 

h n. 

{p)\x-,t)  {py(x,,t) 


(4) 


The  pressure  P is  given  by  a state  equation  {p){x.,t)  = p{{p{x,,t))  and  is  an  artificial  viscosity  proposed  by 
Monaghan  and  Gingold  [9J. 

Two  general  remarks  can  be  add  to  the  presentation  of  the  SPH  method  : 

■ The  continuity  equation  avoids  the  boundary  effects  because  the  density  evolves  only  when  the  distances 
between  particles  change. 

■ Corrections  have  to  be  introduced  in  this  formulation  to  take  into  account  the  space  and  the  time  variations 
of  the  smoothing  length.  The  time  variations  of  h are  described  by  the  law  of  Benz  [lOJ  : 

T'  y(v(x,,t)-v(x  t))-y,W(U,  -xih) 
dt  3 M P\Xj,t)  ' ' 


(5) 
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The  space  variations  of  h are  obtained  by  the  scatter  approximation  : h(x,  y)  = + 

This  relation  should  be  reintroduced  in  all  the  previous  mechanical  equations. 


2.2.  Application  of  the  SPH  method  to  the  bird  strike 

The  projectile  is  modeled  by  a set  of  Lagrangian  particles  arranged  in  a compact  hexagonal  system.  Each 
particles  is  initially  connected  to  twelve  other  particles  but  the  connectivity  is  time  dependant. 


Figure  1 : projectile  mesh 


The  mass  of  each  particle  is  concentrated  at  the  centre  of  the  outer  sphere.  The  bird  material  model  is  based  on 
the  studies  of  J.  Willbeck  [lOJ.  He  has  developed  a synthetic  bird  model  which  is  composed  of  commercial 
gelatin  and  10%  of  phenolic  micro-balloons. 

The  state  equation  is  based  on  the  pressure-density  relationship  of  each  constituent.  But  to  improve  the 
calculation,  this  law  distinguishes  a shock  and  a reduction  in  pressure.  This  model  gives  very  accurate  loading 
including  the  Hugoniot  peak  and  the  stagnation  pressure. 

The  SPH  method  used  the  cubic  kernel  with  compact  support  [1],  [1 IJ  : 


h 


The  contact  algorithm  realize  the  coupling  between  SPH  and  the  classical  Finite  Elements  Methods.  This 
interaction  is  treated  as  the  shock  of  two  spheres  : one  is  a particle  (bird)  and  the  other  is  a factitious  node 
representing  the  nodes  of  an  element  (blade). 


Figure  2 : contact  detection 

The  global  momentum  is  conserved  in  the  normal  direction  whereas  possible  friction  forces  are  not  taken  into 
account  in  the  tangent  direction  [12J  : 

(7) 


e particle 

d defending  node 

tilde  stands  for  ‘before  impact’ 

The  bird  presents  an  hydrodynamics  behaviour  during  impact,  so  the  normal  velocities  after  the  shock  are 
chosen  equal  (shock  without  rebound) : 


v.=v,  (8) 
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These  equations  (7),  (8)  allow  to  calculate  the  velocity  and  the  acceleration  of  the  bird  particles.  This  method 
gives  accurate  results  only  if  the  diameter  of  the  particles  and  the  size  of  the  elements  are  equivalent. 

3.  INDUSTRIAL  APPLICATION 

Bird  strikes  are  today  currently  simulated  on  civil  fan  using  SPH  method.  The  results  for  all  kinds  of  blades 
(wide  chord  solid  blades,  hollow  fan  blades,  shrouded  blades  ...)  are  very  close  to  experimental  data 
[2J,[4J,[12]. 

This  method  has  been  applied  to  anticipate  the  mechanical  behaviour  of  blades  under  different  impact 
conditions  than  fan  blades  : the  blades  dimension  and  the  rotation  speed  are  different  from  fan  civil 
applications.  The  nominal  rotation  speed  is  around  13.600  rpm  and  the  seagull  is  fired  with  an  axial  velocity 
of  197  m/s. 

To  limit  CPU  time,  only  the  impacted  part  of  the  rotor  is  represented.  Each  blade  is  meshed  with  shell 
elements  and  the  disk  with  brick  elements.  The  structural  mesh  approximately  includes  33.000  degrees  of 
freedom  (cf.  figure  3). 

To  compute  accurately  the  bird  slicing,  the  blades  must  be  in  motion  in  the  global  frame.  A procedure  was 
defined  to  describe  the  rotation  of  the  impacted  structure.  This  procedure  follows  three  steps  : 

■ calculation  of  the  centrifugal  stresses 

■ determination  of  the  bird  slicing 

■ deceleration  of  the  rotor. 

The  first  step  is  essential  because  the  dynamic  behaviour  of  the  blades  are  widely  influenced  by  this  additional 
strength. 

The  boundary  conditions  must  take  into  aeeount  the  cyclic  symmetry  of  the  disk  to  model  the  effect  of  the 
rotor  which  is  not  meshed.  The  centrifugal  forces  are  determined  in  the  rotating  frame.  They  are  applied  on 
every  node  and  tend  to  expand  the  structure  in  the  radial  direction.  Thus,  it  is  necessary  to  impose  the 
displacement  continuity  for  each  node  of  the  lateral  disk  faces  located  at  the  same  radial  position  (9). 


radial  velocity 
Vg  tangential  velocity 

(The  mark  1 refers  to  one  face  and  the  mark  2 to  the  other  face) 

The  calculation  is  conducted  until  the  stresses  are  stabilised  (the  structure  is  damped  on  the  first  bending  mode 
of  a blade). 

The  bird  impact  appends  on  the  global  frame  to  represent  correctly  the  slicing  and  the  centrifugal  stresses  of 
the  deformed  structure.  Numerically,  the  rotation  is  introduced  by  prescribed  displacement  on  the  nodes 
located  on  the  lateral  disk  faces.  The  initial  velocities  is  automatically  determined  by  the  code  from  the  radial 
position  of  each  node  and  from  the  rotation  speed  of  the  rotor. 

The  objective  of  this  simulation  is  to  anticipate  the  mechanical  behaviour  of  the  blades,  so  the  outer  shape  of 
the  projectile  does  not  represent  a real  geometry  of  bird  but  a right  circular  cylinder  (cf.  figure  1). 

The  material  behaviour  of  the  blades  is  described  by  a dynamic  elastic-plastic  law.  Contact  conditions  are 
imposed  between  each  blade  and  the  bird  but  also  between  blades  to  avoid  penetrations  under  high  bending 
moments. 

The  results  are  presented  on  a sector  of  three  blades  (figure  3-4).  Snecma  does  not  have  test  data  to  validate 
the  computations,  but  the  first  results  seem  to  predict  realistic  behaviour  of  blades.  The  deformation  under 
impact  is  more  severe  than  impact  on  a civil  engine.  During  the  bird  slice,  larger  deflections  occurred  at  the 
leading  edge.  The  bird  slice  cut  by  each  airfoil  is  not  regular  (even  if  the  bird  density  is  constant)  because 
blades  and  projectile  interact  (figure  4). 

No  active  element  failure  criteria  is  implemented  in  the  code.  The  predicted  maximum  dynamic  strain  and  its 
location  indicate  the  possible  cracking  and  material  loss.  The  value  of  this  criteria  is  determined  with 
experimental  data. 
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This  simulation  provide  useful  data  to  validate  the  design  and  the  load  definitions  under  impact.  The  global 
response  of  the  structure  under  impact  depends  on  the  initial  geometry  and  the  simulation  determines  the 
mode  which  is  excited  (torsion  or  bending).  For  example,  a torsional  deflection  can  lead  to  peak  strain 
response  at  the  airfoil  root  region  and  failure. 

Furthermore,  the  predicted  response  allow  to  validate  the  engine  design.  The  computation  provide  information 
about  the  possible  contact  between  the  impacted  blades  and  : 

■ the  fan  casing  (radial  deflection  - figure  5) 

■ the  following  blades  (tangential  deflection  - figure  6) 

■ the  other  compressor  stage  (axial  deflection  - figure  7) 


4.  CONCLUSIONS 

The  bird  slicing  is  especially  well-predicted  without  any  assumption.  The  results  confirm  the  major  interest  of 
the  SPH  method. 

The  calculations  were  not  validated  by  experimental  data  but  the  global  and  local  deformations  seem  to  be 
well  predicted  : the  blades  interact  with  the  bird  and  the  results  are  coherent  with  the  higher  speed  impact  in 
such  conditions. 

The  SPH  method  has  demonstrated  its  ability  to  be  used  as  an  efficient  analytical  tool  for  predicting  impact 
loads  and  damages  due  to  bird  ingestion.  Since  the  analysis  does  not  require  important  computer  resource  and 
CPU  time,  it  can  be  used  during  the  conception  iterations  and  the  qualification  tests  are  performed  on  the  final 
optimised  geometry.  The  computations  contribute  to  improve  the  quality  of  prediction  of  the  dynamic 
behaviour  of  the  blades  under  impact  and  therefore  to  increase  the  safety  of  flights. 

ACKNOWLEDGEMENTS 

The  authors  would  acknowledge  the  help  of  M.  Lepareux,  M.  Bung  and  M.  Letellier  for  the  implementation  of 
the  SPH  method  in  the  EUROPLEXUS  code. 

REFRENCES 

[IJ  A.  Letellier.  Contribution  a la  modelisation  d’oiseaux  sur  les  aubes  des  reacteurs.  PhD,  Universite 
d’Evry  1996 

[2J  A.  Letellier,  M.  Berthillier,  H.  Bung  and  P.  Galon.  Bird  impact  on  fan  blade  analysis  using  SPH 
method  coupled  with  finite  elements.  Symposium  on  structures  under  extreme  loading  conditions. 
Orlando,  1997 

[3J  A.  Legay,  D.  Leray.  PLEXUS,  introduction  de  la  viscosite  SPH.  Rapport  CEA, 
SEMT/DYN/RT/98.040/A,  1998 

[4J  P.  Vignolles,  P.X.  Bussonet  and  J.  Talbotec.  Conception  des  aubes  fan  soumises  a Timpact  d’oiseaux. 
AGARD  CP  1995 

[5]  W.  Benz.  Numerical  modelling  of  non-linear  stellar  pulsation  : problems  and  prospects.  Kluwer 
Academic,  Boston,  MA,  269-288,  1990 

[6J  G.R.  Johnson,  R.  Stryk.  Dynamic  three  dimensional  computations  for  solids  with  variable  nodal 
connectivity.  Int  J for  Num  Methods  (28),  1989 

[7J  G.R.  Johnson,  E.  Peteren,  R.  Strik.  Incorporation  of  an  SPH  option  into  the  EPIC  code  for  a wide 
range  of  high  velocity  impact  computations.  Int  J of  Impact  Engineering  (14),  1993 
[8J  S.W.  Attaway,  M.W.  Heinstein  and  J.W.  Swegle.  Coupling  of  smooth  particle  hydrodynamics  with 
the  finite  element  method.  Nuclear  Engineering  and  design  (150),  199-205,  1994 
[9J  J.J.  Monaghan,  R.A.  Gingold.  Smoothed  particle  hydrodynamics  : theory  and  applications  to  non- 
spherical  stars,  Montly  Notices  of  the  Royal  Astronomical  Society  (181),  1997 
[lOJ  W.  Benz.  Smooth  Particle  Hydrodynamics  : a review.  Harvard-Smithsoniam  Center  for  Astrophysics, 
1989 

[1 1 J J.  Wilbeck.  Impact  behaviour  of  low  strength  projectiles.  Technical  Report  77-134,  Air  Force  Material 
Laboratory,  Wrigth-Patterson  Air  Force  Base,  Ohio  45433,  1978 
[12J  J.  Bonini,  H.  Bung.  Structural  Impact-Friction  modelling  in  Explicit  finite  element  code  with  the 
Lagrange  Multiplier,  14*  Conference  on  Structural  Mechanics  in  Reactor  Technology,  vol.  6,  37-44, 
1997. 


31-8 


Paper  #31 

Discusser's  Name:  Azzeddine  Soulaimani 
Author's  Name:  D.  Chevrolet 

Q:  How  do  you  apply  the  solid  boundary  conditions  for  the  SPH  method? 
A:  We  are  simply  using  conservation  of  the  momentum  at  the  boundary. 


